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CONTRACTOR REPORT 


THE NASA MARSHALL ENGINEERING 
THERMOSPHERE MODEL 


INTRODUCTION 

The model currently used at NASA/Marshall Space Flight Center to 
describe the properties of the neutral atmosphere between 90 and 
2500 km altitude is the MSFC/J70 model Tl/2], as described in 
[3]. In the Earth Science and Applications Division of the 
Structures and Dynamics Laboratory the computer program used to 
output data from the MSFC/J70 model is the J70MM. This program, 
which runs on an HP-1000 computer, outputs data in a 
latitude/longitude matrix grid (for a fixed time, altitude and 
solar and geomagnetic indices) which can then be printed on a 
line printer. A listing of this program, which is written in 
Fortran IV, can be found in [3]. 

For many purposes it is more convenient to have any computer gen- 
erated output plotted rather than simply tabulated, and with 
several different input parameters used in the model (solar and 
geomagnetic indices, latitude, longitude, altitude and time) it 
was decided that the computer code should be substantially 
modified and streamlined. This was implemented on a VAX 11/780 
computer by the use of more subroutines and by writing the code 
in Fortran 77. Some other modifications were also performed 
which relate to the calculation of ephemeris data and to the cal- 
culation of helium number densities in the upper thermosphere. 
Some additional atmospheric thermodynamic parameters are now also 
calculated in the model. 

In the course of producing global contour plots from the data 
produced by the new streamlined model, it was found that oc- 
casionally density discontinuities were present over certain 
regions of the globe. On further investigation, it was dis- 
covered that these discontinuities were associated with an oc- 
casional problem in the integration routine used in the model. 
More details concerning the nature of this problem and its solu- 
tion are given in this report. 

The final version of this model has been named the NASA Marshall 
Engineering Thermosphere (MET) Model. It consists of a number of 
subroutines which are driven by a very simple driving program, as 
given in the appendix. To produce plots, all one needs to do is 
modify the driving program. This has been done in an accompany- 
ing report [4], the details of which will not appear here. 



SUBROUTINE STRUCTURING 


A schematic representation of the MET model is given in Figure 1. 
The driving program calls the control subroutine ATMOSPHERES, 
which then makes decisions based upon the values of the four 
switches which have been set in the driver. In the basic model 
presented here, only one calculation is performed (ie. calcula- 
tions are performed for one altitude, latitude, longitude and 
time and for one set of solar and geomagnetic conditions) so that 
all of the switches have been set to a character value of "Y" . 
The development of a driving program suitable for the generation 
of data files to be used for creating plots would require the ex- 
ecution of multiple calculations, and hence would require one of 
the switches being set to a character value of "N" in the driver 
program for an x-y line plot, or would require at least one or 
two of the switches being set to a character value of "N" in the 
driver program for a contour plot. The subject of multiple ex- 
ecutions is not discussed any further here, but is discussed in 
[4] . 

Once the subroutine ATMOSPHERES has been called, it successively 
calls the subroutines TIME, SOLSET, J70 and J70SUP. Here, in 
this most basic setup, TIME and SOLSET are called once only. 
Each of these subroutines prompt the user for the time (year, 
month, day, hour and minute) and the solar and geomagnetic in- 
dices (AP, F10.7 and F10.7), respectively, as discussed in [3]. 
In the MET model the user, not given a choice of using the KP in- 
dex, could easily edit SOLSET should the KP index be preferred 
over the use of the AP index. 

Subroutine J70 calculates most of the atmospheric parameters 
which the MSFC/J70 model did, but it contains several modifica- 
tions which will be discussed shortly. J70 successively calls 
the subroutines TME, TINF, JAC, SLV and SLVH. Apart from chang- 
ing the code from Fortran IV to Fortran 77, TINF, SLV and SLVH 
are basically very similar to their old versions appearing in the 
MSFC/J70 model. Some constants have been changed in TME, while 
JAC has been extensively re-written in an attempt to make it in- 
telligible to the reader! An intermittent problem within JAC has 
also been corrected, and this will be discussed in more detail 
shortly. Between altitudes of 440 and 500 km J70 also calls a 
new subroutine called FAIR5, which we will discuss shortly. 

Once calculations in J70 are finished, J70SUP is called. This 
subroutine calculates additional atmospheric thermodynamic 
parameters which will also be described shortly. 
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Figure 1 . The MET MODEL* 


FAIR5 


♦Sequential numerical operation is down and to the 
right 
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SUBROUTINE TME 


A number of relatively minor changes have been made to some of 
the ephemeris constants, as described in [5]. These changes, 
listed in Table 1, lead to improvements in the calculated values 
of the solar declination and solar hour angles, as given in [5]. 
However, the associated changes in the total density are typi- 
cally fractions of a percent. 


TABLE 1 


PARAMETER 

OLD VALUE 

NEW VALUE 

A2 

36000.76854 

36000.76892 

A4 

0.25068447 

0.250684477 

B1 

0.017203 

0.0172028 

B3 

1.410 

1.407 

FMJD 

XMJD-2435839 

XMJD— 2435839+GMT/1440. 

XJ 

(XMJD-2415020 . 0)/36525.0 

(XMJD-2415020. 5)/36525.0 

XLS 

AMOD (Y1+B2*SIN (YI ) -B3 , TPI ) 

AMOD (YI+B2*SIN(Y2)-B3, TPI) 
Where Y2=0 . 017202* (FMJD-3 . ) 

B4 

(TPI/360 . ) *23 . 45 

(TPI/360. ) *(23 .4523-0. 013*XJ) 


SUBROUTINE JAC 

This subroutine has been extensively streamlined and improved by 
the introduction of a function statement (GRAVITY) , two function 
subprograms (TEMP and MOL_WT) and a subroutine (INTEGRATE) . Sec- 
tion 1 of JAC is used for altitudes between 90 and 105 km, while 
Section 2 is used only for altitudes greater than 105 km. Sec- 
tion 3 is used only for the calculation of hydrogen number den- 
sities which below 500 km altitude are set equal to 10 6 m" . 

The function subprogram M0L_WT only calculates the molecular 
weight at and below altitudes of 105 km using the sixth-degree 
polynomial given in [1]. Given this molecular weight and the 
temperature, the total mass density is computed by integrating 
the barometric equation. For altitudes above 105 km the number 
density of each individual specie is calculated by integrating 
the diffusion equation upwards from 105 km. Thus, above 105 km 
altitude the molecular weight is not known (it is calculated 
later from the number densities of all of the species) nor is it 
needed. Thus, the subprogram M0L_WT simply sets the molecular 
weight to unity above 105 km altitude. 
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SUBROUTINE INTEGRATE 


Between 90 and 105 km altitude the atmosphere is well mixed, but 
3S a result of the dissociation of O 2 to 0 the mean molecular 
mass in the model varies in such a manner that the ratio 
n (o)/n(0 2 )=1.5 at 120 km [1]. Above 105 km each individual 
specie is in a state of diffusive equilibrium. The only excep- 
tion to this is hydrogen, which in the model exists above 500 km 
in a state of diffusive equilibrium. Thus, between 90 and 105 km 
the density is computed by integrating the barometric equation 
while above 105 km it is computed by integrating the diffusion 
equation. This integration is performed in the subroutine IN- 
TEGRATE using Simpson's Rule. 


Integration is performed between the altitudes A and D where at 
these altitudes the integrand has values of FA and FD, respec- 
tively. At and below 105 km altitude gM/T is integrated (g is 
the gravitational acceleration, M is the mean molecular weight 
and T is temperature) , while above 105 km altitude g/T is in- 
tegrated. Integration is performed using Simpson's Rule quadra- 
ture, and the value of the integral of either gM/T or g/T between 
altitudes A and D is set equal to R(N) , where the number of sub- 
intervals used to evaluate the integral in the interval A to D is 
given by 2 W . Integration is said to have converged when the 
solution obtained by doubling the number of subintervals is 
"equal" to the previous solution within some specified tolerance, 
e. Mathematically, this is stated as: 


e|R(N)| > | R(N) -R(N-l) | 


( 1 ) 


If equation (1) is satisfied, then convergence is said to have 
been achieved, the result (RR) is set to R(N) , and execution 
returns to the calling program (in this case subroutine JAC) . 
This method is the same as that used in the original J70MM 
program, but the code has been substantially re-written. 


While there is nothing wrong with this basic method, the im- 
plementation of it in J70MM contained a major flaw which appears 
to be intermittent in nature. Figure 2(a) is a contour plot of 
total mass density over a portion of the globe at 230 km al- 
titude. The input conditions required to generate this plot 
are: - 
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Time: 

Indices: 


21 June, 1987, 0400 hrs. 
f 10.7 = 230 ' f 10.7 = 230 ' A p = 0 

Altitude: 230 km 


This plot was not generated with the MET model, but with a 
modified form of the original J70MM program. The important fea- 
ture of this plot is not the actual magnitudes of the densities, 
but the discontinuous nature of these densities as a function of 
both latitude and longitude, as revealed by the several "saw- 
toothed curves" which are clearly evident in the plot. It is im- 
portant to note, however, that all of the empirical equations 
upon which the model is based are continuous functions of both 
latitude and longitude [1]. Also, each individual data point 
used in Figure 2(a) is independent of all of the other data 
points which are used. Hence, any error in the model which has 
produced the discontinuities in this figure are peculiar to the 
calculations made for each individual data point which is in er- 
ror. Due to the fact that a contour plot of mass density at 200 
km altitude, (with all other input parameters identical to those 
used in Figure 2(a)) appeared to be perfectly continuous over 
latitude and longitude, a discontinuity in density as a function 
of altitude was sought. 

Figure 2(b) shows mass density as a function of altitude for the 
same time and solar and geomagnetic indices as previously used in 
Figure 2(a), and for a latitude of -62.1° and a longitude of 
295°. [This geographic position was chosen because it cor- 
responds to a region of discontinuity in Figure 2(a). Also, only 
a small altitude range centered about 230 km was considered. Had 
a large altitude range been chosen, any small discontinuities may 
have been hidden on the plot by the scaling required to display 
the rapid exponential decrease of density which occurs over in- 
creasing altitude. ] Immediately evident in this figure is a den- 
sity discontinuity centered at approximately 229.8 km altitude 
with a vertical extent of approximately 1 km and a magnitude of 
about 1.2%. This discontinuity is a direct result of the conver- 
gence criterion used in the integration procedure, and although a 
detailed description of this problem is beyond the scope of this 
report, a brief description of the problem and its solution fol- 
lows. 

As already stated, integration is performed using Simpson's Rule 
quadrature. The form of this procedure, as used here, is an 
iterative one, and is such that unnecessary re-evaluation of the 
integrand is not required at the "old" subinterval end- and mid- 
points when the number of subintervals is doubled. Specifically, 
if the number of subintervals is 2 , then the area A(J) is given 
by 
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A ( J) = 1/3 (SONE)j + 4/3hj (STOW)j 


( 2 ) 



ALTITUDE ( KM ) LATITUDE 



LONGITUDE 

Figure 2(a) Contour plot of J70MM model mass density 



- 10 - 3 

DENSITY ( x 10 kp.m ) 

Figure 2(b) Altitude plot of J70MM model mass density 




where hj is the subinterval length (equal to (D-A).2 J in sub- 
routine INTEGRATE) and where SONE and STOW are functions of the 
integrand evaluated at subinterval end- and mid-points (the exact 
functional forms of SONE and STOW are not required here, but can 
be found in subroutine INTEGRATE) . If the number of subintervals 
is now doubled (ie. 2 J+1 ) , then the area calculated with this new 
number of subintervals (A(J+1) ) in terms of A(J) is 


A ( J+l) = 1/2 A ( J) - 2/3 h J+1 (STOW)j + 4/3 h J+1 (STOW) J+1 (3) 


Evaluation of A(J+1) simply requires the evaluation of (STOW) J+1 , 
because A(J) and (STOW)j are already known. Evaluation of 
(STOW) j + -, requires calculating the integrand at every point lying 
halfway "between every "old” subinterval (2 J of them) end- and 
mid-point. 

Examination of equation (3) reveals that if 


4/3 h J+1 (STOW) J+1 - 2/3 h J+1 (STOW) j % 1/2 A(J) (4) 

then iterative convergence, as defined in equation (1) , may be 
achieved if the equality in equation (4) is "sufficiently close" 
(this will, of course, depend on the value of e ) . That equation 
(4) , and hence that equation (1) , could be satisfied by fluke is 
not beyond the realm of possibility, and in fact as it 
transpires, the discontinuities in Figures 2(a) and 2(b) are 
directly attributable to this fluke convergence. For altitudes 
below about 229.3 km and above about 230.3 km the integration 
procedure required that at least 2 5 (ie. 32) subintervals be used 
to ensure convergence. Between these two altitudes only 2 3 (ie. 
8) subintervals were required for equation (1) to be satisfied, 
meaning that no significant increase in accuracy was achieved by 
increasing the number of subintervals from 4 to 8. However, a 
re-calculation using 16 subintervals showed that convergence had 
not been achieved. In fact, 32 subintervals were required to en- 
sure convergence. 

The remedy to this problem, which was implemented in, and subse- 
quently incorporated into INTEGRATE, is a very simple one but may 
not be guaranteed under all conditions. Quite simply, above 180 
km altitude, the number of subintervals used must be at least 16 
(ie. 2 4 ) before convergence is tested using equation (1) . Up un- 
til the present time this has always been successful, but any fu- 
ture encounters with problems related to this should be reported 
to either the author or to the Branch Chief, ED44, NASA/MSFC, 
Huntsville, Alabama 35812. 

The results of using the new MET model, which incorporates these 
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modifications, are shown in figures 3(a) and 3(b). The input 
parameters used to generate these results are identical to those 
used for the corresponding results (using a modified J70MM 
program) shown in Figure 2. The MET model shows no sign of any 
discontinuities . 


SUBROUTINE FAIR5 

In the MET model, as in the J70MM model, the helium number den- 
sity (and concomitantly the total mass density) above 500 km al- 
titude is modified by seasonal-latitudinal variations according 
to the equations of Jacchia [2]. This seasonal-latitudinal 
variation of helium is calculated in subroutine SLVH. (This 
should not be confused with the seasonal-latitudinal variation of 
density which is calculated for altitudes between 90 and 170 km 
in subroutine SLV) . In the J70MM model the inclusion of 
seasonal- latitudinal variations at altitudes above 500 km led to 
a discontinuity in the helium number densities (and also, there- 
fore, in the total mass density) across the 500 km level. This 
undesirable feature of the J70MM model has been eliminated in the 
MET model by the use of subroutine FAIR5 [7]. FAIR5 is invoked 
between 440 and 500 km altitude, while SLVH is now invoked at all 
altitudes above 440 km (instead of 500 km) . Using a spline 
routine, FAIRS adds progressively more and more of the seasonal- 
latitudial variation to the helium number densities as one 
progresses from 440 to 500 km altitude, so that no effect of SLVH 
is included at 440 km altitude, but the full affect of SLVH is 
included at 500 km altitude. While FAIR5 now corrects the den- 
sity discontinuity problem at 500 km altitude, it does not in any 
way significantly effect the total mass density between the al- 
titudes of 440 and 500 km, because helium is still a minor 
species in this altitude range. 


SUBROUTINE J70SUP. 

This subroutine calculates supplementary atmospheric parameters 
which are not directly calculated in J70. The input parameters 
for J70SUP are altitude and the entire OUTDATA array from J70. 
The output parameters calculated by J70SUP are g (gravitational 
acceleration) , H (pressure scale-height) , Cp (specific heat at 
constant pressure), C v (specific heat at constant volume) and 7 
(ratio Cp/C v ) . 

For consistency g is calculated as in [1], given by 

g = 9.80665 (1 + Z/R e ) ~ 2 ms" 2 (5) 


where Z is altitude (in km) and the Earth's radius Rg = 6.356766 
x 10^ km. The equation for the pressure scale-height, H, can be 
found in most standard texts (eg. [6]) and in a hydrostatic and 
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ALTITUDE 







isothermal atmosphere is given by 


H = p/pg (6) 


where p is pressure and p is mass density. The value of 7 is 
then calculated from the weighted mean of all diatomic and 
monatomic gases by assuming that all diatomic gases have 7 = 1.4 
and all monatomic gases have 7 = 1.67. By using equation (6) , 
the equation of state and the fact that R = c_-c v (R is the gas 
constant) one obtains p 


C v = gH/ (7-1) T 


(7) 


where T is temperature. 


TEST CASE 

The following is a test case for the MET model. For the follow- 
ing input: 

Time: 1987, October 27, 1400 hrs 

Indices: F 10>7 = 230, F 1(K7 = 230, AP = 400 

Altitude = 200 km, Latitude = Longitude = 0° 


The following output should obtain: 


Exospheric temperature 

= 

2105.903 

K 


Local temperature 

= 

1495.590 

K 

10 15 

N 2 number density 

= 

6.378800 

X 

0 2 number density 

= 

7.158251 

X 

10 14 

0 number density 

= 

4.897181 

X 

10 lb 

A number density 

= 

1.146983 

X 

io* 3 

He number density 

= 

9.737494 

X 

10 12 

H number density 

= 

1.00000 

X 

10 6 

Average molecular weight 
Total mass density 

= 

23.34522 

4.656589 

X 

10~ 10 

Log 10 (mass density) 
Total pressure 

= 

-12.33193 


= 

2.480329 

X 

10 4 


-3 

-3 

-3 

-3 

-3 

-3 


kgm 


Pa 


-3 


Local gravitational acceleration 
Ratio of specific heats 
Pressure scale-height 
Specific heat at constant pressure 
Specific heat at constant volume 


9.217513 

1.510544 

57786.66 

1053.729 m^s_ 2 K" 1 
697.5829 m 2 s 2 K 1 


For the same input conditions the original J70MM model produced 
an exospheric temperature of 2105. 775K, a local temperature of 
1495. 539K and a mass density of 4.656528 x 10 -10 kgm -3 (the en- 
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tire output from J70MM is not given here) . Therefore, for this 
particular case the differences between results from the MET 
model and the J70MM model are typically less than 0.01%. One 
must remember, however, that occasionally results from the two 
models will differ by more than 1% (eg. Figure 2) . 


CONCLUSIONS AND RECOMMENDATIONS 

This report has described how the MSFC/JO model was modified in 
order to produce the NASA/MET model. Although this report has 
been primarily concerned with the several modifications which 
directly affect the numerical output of the MET model, most of 
the modifications to the MSFC/J70 model entailed extensive 
rewriting of the FORTRAN code with the result that the MET model, 
written is FORTRAN 77, is easy to follow and easy to tailor to 
suit ones individual needs. The MET model computer program, 
listed in the Appendix, should be compared with the J70MM model 
computer program, which is listed in [3], 

Most of the numerical differences between the output of the two 
models is generally small, and typically fractions of a percent. 
Occasionally, however, differences can be about one percent (in 
total mass density) , and these are associated with the intermit- 
tent integration error in the MSFC/J70 model. The MET model out- 
puts several parameters which the MSFC/J70 (or J70MM) model did 
not. These extra parameters may be useful to some users. 

The MET model is currently available to users who have access to 
the NASA/GSFC ENVIRONET network. In the present report, it is 
recommended that the MET model replace the MSFC/J70 and J70MM 
models for all studies relating to all design, building and 
operations of NASA space vehicles. 
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APPENDIX 


C*' 

c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c» 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c» 
c* 
c* 
c* 
c* 
c* 
c* 
c* 

REAL*4 INDATA (12) . OUTDATA (12) . AUXDATA (5) 

CHARACTER *1 SWITCH (4) 

CALL LIB$INIT_TIMER 

C Set all switches to ’Y’ so that only one particular calculation is performed 

SWITCH (1) - ’Y’ 

SWITCH (2) - *Y* 

SWITCH (3) - *Y* 

SWITCH (4) - 'Y* 

CALL ATMOSPHERES ( INDATA, OUTDATA. AUXDATA, SWITCH ) 

C Now type output data 

Type All output in MKS units’ 

Type ’ ’ 

Type »,’ Exospheric temperature *■ ’, OUTDATA (1).’ K* 

Type Temperature » OUTDATA (2),’ K’ 

Type *,’ N2 number density «■ ’, OUTDATA (3),’ /m3’ 

Type 02 number density ” ’, OUTDATA (4),’ /m3’ 

Type *,’ 0 number density • ’, OUTDATA (5),’ /m3’ 

Type *,’ A number density « OUTDATA (6),’ /m3’ 

Type He number density « ’, OUTDATA (7),’ /m3’ 

Type H number density — OUTDATA (8),’ /m3’ 

Type Average molecular wt . « ’, OUTDATA (9) 

Type *,’ Total mass density - OUTDATA (10),’ kg/m3’ 

Type *,’ Log10 mass density - ’, OUTDATA (1l) 

Type Total pressure ” ’, OUTDATA (12),’ Pa’ 

Type *,’ Local grav. acceln. « ’, AUXDATA (1),’ m/sec’ 

Type *,’ Ratio specific heats - ’, AUXDATA (2) 

Type *,’ Pressure scale-height - ’, AUXDATA (3),’ m’ 

Type *,’ Specific heat cons, p — ’ , AUXDATA (4),’ m2. sec-2. K— 1 ’ 


preceding page blank not filmed 


The Marshall Space Flight Center 
Marshall Engineering Thermosphere Model 


written by 
Mike Hickey 

Universities Space Research Association 
NASA / MSFC , ED44 
Tel. (205) 544-5692 

This program is a driving program for the following subroutines :— 

ATMOSPHERES 

SOLSET 

TIME 

J70 

The atmospheric model is a modified Jacchia 1970 model and is given in 
the subroutine J70. All of the other subroutines were designed to 
allow flexible use of this model so that various input parameters could 
be varied within a driving program with very little software development 
Thus, for example, driving routines can be written quite easily to 
facilitate the plotting of output as line or contour plots. Control is 
achieved by setting the values of four switches in the driving program, 
as described in subroutine ATMOSPHERES. 




SUBROUTINE ATMOSPHERES ( INDATA, OUTDATA, AUXDATA, SWITCH ) 


C*< 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 
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c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 


DESCRIPTION:- 


Colculate atmospheric data in single precision using subroutine J70 
and J70SUP. 


SUBROUT I NES:- 


TIME, SOLSET, GMC, J70 and J70SUP 


INPUT :- 


all single precision, either through 
subroutines or from main driver prog. 


INDATA (1) — altitude - Z 

(2) — latitude - XLAT 

(3) — longitude - XLNG 

. . (4) — year (yy) - IYR 

(5) — month (mm) ■ MN 

(6) — day (dd) - IDA 

(7) — hour (hh) - I HR 

(8) — mins (mm) - MIN 

.. (9) — geomagnetic index — IGEO_IND 

.. (10) — solar radio noise flux- F10 

(11) — 162-day average F10 — F10B 


(12) — geomagnetic activity index - GI-AP 


OUTPUT :- 


NOTE : All output in MKS units 


all single precision 


OUTDATA 


(1) — exospheric temperature (K) 

(2) — temperature at altitude Z 

(3) — N2 number density (per meter-cubed) 


(4) — 02 number density ( .. ) 

(5) — 0 number density ( .. ) 

(6) — A number density ( .. ) 

(7) — He number density ( .. ) 

(8) — H number density ( .. ) 


(9) — average molecular weight 

(10) — total density 

(11) — Iog10 ( total density ) 

(12) — total pressure ( Pa ) 


AUXDATA (1) — 
(2) — 

(3) — 

(4) — 

(5) — 


gravi tat ional acceleration ( m/s-s ) 
ratio of specific heats 
pressure scale— height ( m ) 
specific heat at constant pressure 
specific heat at constant volume 


COMMENTS :- 


SWITCH(I) — if Y(es) , 
SWITCH(2) — if Y(es) , 
SWITCH(3) — if Y(es) , 
SWITCH(4) — if Y(es) , 


date and time are input from terminal through 
subroutine TIME once only 

solar/magnetic activity are input from terminal 
through subroutine SOLSET once only 
only ONE altitude value is input from terminal 
through main calling program 

only (WE latitude AND longitude are input from 
terminal through main calling program 


C* ATMOSPHERES written by Mike Hickey ( USRA, NASA/ED44 ) • 

C* Tel: (205) 544-5692 • 

C* January-Apr i I 1987 * 

C***»* * ** ****** ****************** 


EXTERNAL TIME 
DIMENSION AUXDATA (5) 

INTEGER HR 

REAL*4 LAT, LON. INDATA (12). OUTDATA (12) 
CHARACTER *1 SWITCH (4) 

PARAMETER PI - 3.14159265 


C — — — 

C Thie next section is only executed on the first call to ATMOSPHERES 

DO WHILE ( CALL. EQ. 0.0 ) 

C SECTION A:- 
C 

IF ( SWITCH(I). EQ. *Y* ) THEN 

CALL TIME ( IYR, MON, IDA, HR. MIN. SWITCH(I) ) 

INDATA (4) - FLOAT J (IYR) 

INDATA (5) - FLOAT J (MON) 

INDATA (6) - FLOAT J (IDA) 

INDATA (7) - FLOAT J (HR) 

INDATA (8) - FLOAT J (MIN) 

END IF 

C SECTION B:- 
C 

IF ( SWITCH(2) . EQ. *Y* ) THEN 

CALL SOLSET ( IGEO_IND, F10, F10B, GI. SWITCH(2) ) 

INDATA (9) - FLOAT J (IGEO_IND) 

INDATA (10) - F10 
INDATA (11) - F10B 
INDATA (12) - GI 

EM) IF 

C SECTION C:- 
C 

IF ( SWITCH(3) . EQ. *Y' ) THEN 

TYPE *.* Input altitude, km’ 

ACCEPT • , INDATA (1) 

Z - INDATA (1) 

ENO IF 

C SECTION D:- 
C 


IF ( SWITCH(4) . EQ. ’Y* ) THEN 

TYPE *,* Input latitude and longitude, degrees' 

ACCEPT *. ( INDATA(I) , I- 2,3 ) 

LAT - INDATA (2) 

LON - INDATA (3) 

RLT - INDATA (2) • PI / 180. I geographic latitude, radians 
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END IF 



o o 


CALL - 1.0 

END DO 

End of first executable section 


C The following depend on the values of the switches 


C**«* 

C* SECTION 1:- 


IF ( SWITCH(I). NE. *Y* ) THEN 


IYR - JINT ( 
MON - JINT ( 
IDA - JINT ( 
HR - JINT ( 
MIN - JINT ( 
CALL TIME ( 


INDATA (4) ) 
INDATA (5) ) 
INDATA (6) ) 
INDATA (7) ) 
INDATA (8) ) 

IYR. MON. IDA. HR, 


MIN. SWITCH(I) ) 


END IF 


C* SECTION 2:- 


IF ( SWITCH(2) . NE. ’Y - ) THEN 

IGEO_IND - JINT ( INDATA (9) ) 

F10 - INDATA (10) 

F10B - INDATA (11) 

GI - INDATA (12) 

CALL SOLSET ( IGEO_IND, F10. F10B. GI. SWITCH(2) ) 
END IF 


C* SECTION 3:- 


IF ( SMTCH(3). NE. *Y* ) THEN 
Z - INDATA (1) 

END IF 


C***«* 

C* SECTION 4:- 


IF ( SWITCH(4). NE. *Y' ) THEN 

LAT - INDATA (2) 

LON - INDATA (3) 

RLT - INDATA (2) • PI / 180. I geographic latitude, rad 
END IF 


C All setting-up complete. 

CALL J70 ( INDATA, OUTDATA ) 

CALL J70SUP ( Z. OUTDATA, AUXDATA ) 


RETURN 

ENTRY ATMOS_ENT ( DUMMY ) 

CALL - DUMMY 

RETURN 

END 


i ans 
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o o o 


SUBROUTINE TIME ( IYR, MON, IDA. HR, MIN, SWITCH ) 


c* 



c* 


This subroutine sets up time 

c* 



c* 


INPUTS/OUTPUTS: 

c* 



c* 

IYR 

« year ( 2 digits ) 

c* 

MON 

- month 

c* 

IDA 

- day of month 

c* 

HR 

» hour of day 

c* 

MIN 

*■ minutes 

c* 



c* 


Written by Mike Hickey, USRA 


DIMENSION IDAY ( 12 ) 

INTEGER HR 
CHARACTER* 1 SWITCH 

DATA IDAY/ 31. 28, 31, 30, 31. 30, 31. 31. 30, 31, 30. 31 / 
PARAMETER PI - 3.14159265 


If SWITCH - Y(es) then input data and time from terminal 


IF ( SWITCH. EQ.’Y’. OR. SWITCH. EQ. ’y’ ) THEN 

TYPE «, ' Input date and time of date? ( yy ,mm,dd,hh,mm ) ’ 
ACCEPT *, IYR, MON, IDA. HR, MIN 

END IF 


IF ( JMOD (IYR. 4) .EQ. 0 ) THEN 

IF ( JMOD (IYR, 100) .NE. 0 ) IDAY ( 2 ) - 29 

ELSE 

IDAY ( 2 ) - 28 

END IF 


DAYTOT - 0.0 


DO 1 I - 1 , 12 

DAYTOT - DAYTOT + FLOAT J ( IDAY ( I ) ) 

1 CONTINUE 


IF ( MON. GT. 1 ) THEN 

KE - MON - 1 
ID - 0 

DO 2 I - 1, KE 
ID - ID + IDAY (I) 

2 CONTINUE 


ID - ID + IDA 
DD - IDA 

END IF 


ELSE 
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RETURN 

END 



SUBROUTINE SOLSET ( IGEO.IND, F10, F10B, GI, SWITCH ) 


C****** ******* •••***»*•*** •*** ***** 

c* 

C* This subroutine simply calls for a setup of the solar-act ivty and auroral 
C* activity indices. 


c* 

c* 


INPUTS/OUTPUTS: 

c* 

c* 

IGEO_IND - 

geomagnetic index 

c* 

FI 0 

solar radio noise flux 

c* 

F10B 

162— day average F10 

c* 

GI 

geomagnetic activity index 

c* 

c* 

Wri tten 

by Mike Hickey, USRA 

c** 




CHARACTER* 1 SWITCH 
IGEO.IND - 2 

C 

C If SWITCH - Y(es) then input geomagnetic indices from terminal 

C 

IF ( SWITCH. EQ.’Y*. OR. SWITCH . EQ. ’ y ’ ) THEN 

C TYPE *, ' Input geomagnetic index ( 1-KP, 2-AP ) * 

C ACCEPT *. IGEO.IND 

TYPE *, * Input solar radio noise flux ( F10 « 0-400 ) ’ 

ACCEPT *. F10 

TYPE *, ' Input 162-day average F10 ( F10B - 0-250 ) ’ 

ACCEPT *, F10B 

C IF ( IGEO.IND . EQ. 2 ) 

C 

C TYPE *, ’ Input geomagnetic activity 

C 
C 
C 

C TYPE *. ’ Input geomagnetic activity 

C 

C END IF 

TYPE *,* Input A P index ( AP - 0 - 400 ) ’ 

ACCEPT *, GI 

END IF 

C 


THEN 

index ( GI - 0-400 ) ’ 
ELSE 

index ( GI « 0-9 ) ’ 


RETURN 

END 
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oooooooooooooooooooooooooooooooooo 


SUBROUTINE J70SUP ( Z, OUTDATA, AUXDATA ) 


DESCRIPTION:- 


J70SUP calculates auxilliary variables which are output in array 
AUXDATA, given data input from J70 which are contained in array OUTDATA. 

INPUT DATA:- 


z — 

altitude (km) 



TZZ — 

temperature at altitude z 

- 

OUTDATA (2) 

— 

N2 number density 


.. (3) 

— 

02 . . 

■ 

•• (4) 

— 

0 

- 

.. (5) 

— 

A 

- 

(6) 

— 

He . . 

- 

.. (7) 

— 

H 

- 

•• (8) 

EM — 

average molecular weight 

• 

• • 0) 

DENS — 

total density 

■ 

.. (10) 

P — 

total pressure 

- 

.. (12) 


OUTPUT DATA:- 


G — gravitational acceleration 
GAM — ratio of specific heats 
H — pressure scale-height 
CP — specific heat at constant pressure 
CV — specific heat at constant volume 


AUXDATA (1) 
AUXDATA (2) 
AUXDATA (3) 
AUXDATA (4) 
AUXDATA (5) 


Written by Mike Hickey. USRA 


REAL*4 OUTDATA (12), AUXDATA (5). H 

G - 9.80665 / ( ( 1. + Z / 6.356766E3 )**2 ) 

H - OUTDATA (12) / ( G * OUTDATA (10) ) 

SUM1 - OUTDATA (3) + OUTDATA (4) 

SUM2 - 0.0 

DO 1 I - 5, 8 
SUM2 - SUM2 + OUTDATA (I) 

1 CONTINUE 

GAM - ( 1.4 • SUM1 + 1.67 • SUM2 ) / ( SUM1 + SUM2 ) 

CV - G • H / ( ( GAM - 1 .0 ) • OUTDATA (2) ) 

CP - GAM • CV 

AUXDATA (1) - G 
AUXDATA (2) - GAM 
AUXDATA (3) - H 
AUXDATA (4) - CP 
AUXDATA (5) - CV 

RETURN 

END 
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o o o 


SUBROUTINE J70 ( INDATA, OUTDATA ) 


c** 

**************************************************************************** 

c** 





*♦ 

c** 



J70 developed from J70MM by 

** 

c** 



Mike P. Hickey 


** 

c*» 



Universities Space Research 

Association 

** 

c** 



at 


** 

c** 



NASA /Marshall Space Flight Center, ED44, 

** 

c** 



Huntsville, Alabama, 35812, USA. 

** 

c** 



Tel. (205) 544-5692 

** 

c** 





** 

c** 

INPUTS: 


through the subroutine calling list 

** 

c** 





** 

c** 

OUTPUTS 


through the subroutine ca 

ling list 

** 

c** 





♦ * 

c** 





** 

c»» 



INPUT DATA: 


** 

c»* 





** 

c** 


z 

— altitude 

- INDATA (1) 

* * 

c** 


XLAT 

— latitude 

- INDATA (2) 

** 

c** 


XLNG 

— longitude 

- INDATA (3) 

** 

c*» 


IYR 

— year (yy) 

- INDATA (4) 

** 

c** 


MN - 

- month (mm) 

- INDATA (5) 

** 

c** 


IDA 

— day (dd) 

- INDATA (6) 

** 

c** 


I HR 

— hour (hh) 

- INDATA (7) 

• * 

c** 


MIN 

— mins (mm) 

- INDATA (8) 

** 

c** 


11 - 

- geomagnetic index 

- INDATA (9) 

** 

c** 


F10 

— solar radio noise flux 

- INDATA (10) 

** 

c*» 


F10B 

— 162-day average F10 

- INDATA (11) 

** 

c** 


GI - 

- geomagnetic activity index - INDATA (12) 

* * 

c** 





*♦ 

c*» 





** 

c** 



OUTPUT DATA: 


** 

c** 





*# 

c** 


T— 

exospheric temperature 

- OUTDATA (1) 

** 

c** 


TZZ— 

- temperature at altitude Z 

- OUTDATA (2) 

** 

c** 


A ( 1 ) 

— N2 number density 

- OUTDATA (3) 

* * 

c** 


A(2) 

— 02 number density 

- OUTDATA (4) 

** 

c** 


A(3) 

— 0 number density 

- OUTDATA (5) 

* * 

c** 


A(4) 

— A number density 

- OUTDATA (6) 

** 

c** 


A(5) 

— He number density 

- OUTDATA (7) 

** 

c** 


A(6) 

— H number density 

- OUTDATA (8) 

** 

c** 


EM— 

average molecular weight 

- OUTDATA (9) 

** 

c** 


DENS 

— total densi ty 

- OUTDATA (10) 

** 

c** 


DL— 

Iog10 ( total density ) 

- OUTDATA (11) 

** 

c*« 


P— 

total pressure 

- OUTDATA (12) 

** 

c** 





*♦ 

c** 

NB. 

Input 

through array ’INDATA’ 


* * 

C*» 

c** 

******** 

Output through array ’OUTDATA* 


** 
► ** 


DIMENSION A ( 6 ) 

REAL*4 INDATA ( 12 ), OUTDATA ( 12 ) 

PARAMETER RGAS - 8.31432E3 1 J/kmol-K 

PARAMETER BFH - 440.0 

C Calcultions performed for only one latitude , one longitude 
C and one altitude 


parameters to 

INDATA 

va 1 ues 

Z 


INDATA 

(D 



XLAT 

- 

INDATA 

(2) 



XLNG 

- 

INDATA 

(3) 



IYR 

- 

JINT ( 

INDATA 

(4) 

) + 1900 

MN 

- 

JINT ( 

INDA'A 

<5) 

) 

IDA 

- 

JINT ( 

INDATA 

(6) 

) 

I HR 

- 

JINT ( 

INDATA 

(7) 

) 

MIN 

m 

JINT ( 

INDATA 

(8) 

) 

11 

- 

JINT ( 

INDATA 

(9) 

) 

F10 

- 

INDATA 

(10) 




F10B - INDATA (11) 
GI - INDATA (12) 


CALL TME ( MN , IDA , IYR . I HR . MIN , XLAT . XLNG . SDA . 

SHA , DD , DY ) 

CALL TINF ( F10 . F10B , GI , XLAT . SDA , SHA , DY , II . TE ) 

CALL JAC ( Z , TE . TZ . A(1) , A(2) . A(3) . A(4) . A(5) . A(6) , 
EM , DENS . DL ) 


DENLG - 0. 

DUMMY - DL 
DEN - DL 

IF (2 .LE. 170. ) THEN 

CALL SLV ( DUMMY . Z . XLAT , DD ) 
DENLG - DUMMY 

END IF 


C 

C** 'Fair* helium number density between base fairing height ( BFH ) and 500 km 
C 


IF ( Z. GT. 500. ) THEN 

CALL SLVH ( DEN . A(5) . XLAT . SDA ) 

DL - DEN 

ELSE IF ( Z .GT. BFH ) THEN 

DHEL1 - A ( 5 ) 

DHEL2 - A ( 5 ) 

DLG1 - DL 
DLG2 - DL 

CALL SLVH ( DLG2 . DHEL2 . XLAT , SDA ) 

IH - Z 

CALL FAIR5 ( DHEL1 . DHEL2 . DLG1 . DLG2 . IH , FDHEL , FDLG ) 
DL - FDLG 
A ( 5 ) - FDHEL 

END I F 


DL - DL + DENLG 

DENS - 10. **DL 

XLAT - XLAT • 57.29577951 


C Fill OUTDATA array 

OUTDATA (1) - TE 
OUTDATA (2) - TZ 

DO 80 I - 1, 6 

OUTDATA (1+2) - 1.E6 * ( 10. •* A(I) ) 
80 CONTINUE 

OUTDATA (9) - EM 
OUTDATA (10) - DENS * 1000. 

OUTDATA (11) - DL 
P - OUTDATA (10) • RGAS • TZ / EM 
OUTDATA (12) - P 

RETURN 

END 


24 



SUBROUTINE TME ( MN . IDA , IYR . I HR , MIN . XLAT . XLNG , 
SDA , SHA . DD , DY ) 


***** ******* ************************ 

C** Subroutine ’TME’ performs the calculations of the solar declination * 

C** angle and solar hour angle. * 

C** * 

C** INPUTS: MN - month * 

C** IDA - day * 

C** IYR - year • 

C** IHR - hour * 

C«* MIN - minute * 

C** XMJD— mean Julian date * 

C** XLAT- latitude ( input-geocentric latitude ) * 

C** XLNG— longitude ( input-geocentric longitude, -180.+180 ) * 

C** * 

C** OUTPUTS: SDA - solar declination angle (rad) * 

C** SHA — solar hour angle (rad) * 

C** DD - day number from 1 JAN. * 

C** DY - DD / tropical year * 

C** Modified by Mike Hickey, USRA * 

C* ****••***•****•**•*•*••*••••*••***••**••*••*•*•**•*********••**•**********•* 

DIMENSION IDAY(12) 

DATA IDAY / 31,28,31 ,30,31,30 ,31,31,30 ,31,30,31 / 

PARAMETER YEAR - 365.2422 

PARAMETER A1 - 99.6909833 , A2 - 36000.76892 

PARAMETER A3 - 0.00038708 , A4 - 0.250684477 

PARAMETER B1 - 0.0172028 , B2 - 0.0335 , B3 - 1.407 

PARAMETER PI - 3.14159265 , TPI - 6.28318531 

PARAMETER PI2 - 1.57079633 , PI32 - 4.71238898 

PARAMETER RAD_DEG - 0.017453293 

XLAT - XLAT / 57.29577951 
YR - IYR 

IF ( JM0D(IYR,4) .EQ. 0 ) THEN 

IF ( JMOD(IYR,100) .NE. 0 ) IDAY(2) - 29 I Century not a leap year 
ELSE 

IDAY(2) - 28 
END IF 

ID - 0 

IF ( MN. GT. 1 ) THEN 
DO 20 I - 1 , MN— 1 
ID - ID + IDAY(I) 

20 CONTINUE 

END IF 

ID - ID + IDA 
DD - ID 
DY - DDAEA R 
C 

C** Compute mean Julian date 
C 

XMJD - 2415020. + 365. * ( YR - 1900. ) + DD 

+ FLOAT J ( (IYR - 1901 ) / 4 ) 

C 

C** Compute Greenwich mean time in minutes GMT 
C 

XHR - IHR 
XMIN - MIN 

GMT - 60 • XHR + XMIN 
FMJD - XMJD - 2435839. + GMT / 1440. 

C 

C** Compute Greenwich mean position - GP ( in rad ) 

C 

XJ - ( XMJD - 2415020.5 ) / ( 36525.0 ) 

GP - AMOD ( A1 + A2 • XJ + A3 • XJ • XJ + A4 • GMT , 360. ) 

C 

C** Compute right ascension point - RAP ( in rad ) 

C 

C** 1st convert geocentric longitude to deg longitude - west neg , + east 
C 
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IF ( XLNG .GT. 180. ) XLNG - XLNG - 360. 
RAP - AMOD ( GP + XLNG . 360. ) 


C 

C** Compute celestial longitude - XLS ( In rad ) zero to 2PI 

C 

Y1 - B1 • FMJD 

Y2 - 0.017202 * ( FMJD - 3. ) 

XLS - AMOD ( Y1 + B2 * SIN(Y2) - B3 . TPI ) 

C 

C** Compute solar declination angle - SDA ( in rad ) 

C 

BA - RAD.DEG • ( 23.4523 - 0.013 * XJ ) 

SDA - ASIN ( SIN ( XLS ) • SIN ( B4 ) ) 

C 

C** Compute right ascension of Sun - RAS ( in rad ) zero to 2PI 

C 

RAS - ASIN ( TAN ( SDA ) / TAN ( B4 ) ) 

C 

C»* Put RAS in same quadrant as XLS 
C 

RAS - ABS ( RAS ) 

TEMP - ABS ( XLS ) 

IF ( TEMP. LE. PI .AND. TEMP.GT.PI2 ) THEN 
RAS - PI - RAS 

ELSE IF ( TEMP.LE.PI32 .AND. TEMP. GT. PI ) THEN 
RAS - PI + RAS 

ELSE IF ( TEMP.GT.PI32 ) THEN 
RAS - TPI - RAS 

END IF 

IF ( XLS. LT. 0. ) RAS - -RAS 
C 

C** Compute solar hour angle - SHA ( in deg ) 

C 

SHA - RAP • RAD.DEG - RAS 

IF ( SHA.GT.PI ) SHA - SHA - TPI 

IF ( SHA. LT .—PI ) SHA - SHA + TPI 


RETURN 

END 
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SUBROUTINE TINF ( F10 , F10B , GI. XLAT, SDA , SHA . DY , II . TE ) 


C******* ****** 

C** Subroutine ’TINF’ calculates the exospheric temperature according to 

C** L. Jacchia SAO 313. 1970 

C** 

C** F10 - solar radio noise flux ( x E-22 Watts / m2 ) 

C** F10B- 162-day average F10 

C** GI - geomagnetic octivity index 

C** LAT - geographic latitude at perigee ( in rad ) 

C** SDA « solar declination angle ( in rad ) 

C** SHA » solar hour angle 

C** DY - D / Y ( day number / tropical year ) ; 1 

C** II - geomagnetic equation index ( 1 — GI-KP , 2 — GI«AP ) 

C** RE « diurnal factor KP, F10B, AVG 

C** 

C** CONSTANTS — C - solar activity variation 

C** — BETA , etc - diurnal variation 

C** — D « geomagnetic variation 

C** — E » semiannual variation 

C** 

C*« Modified by Mike Hickey, USRA 

C**************************************************************************** 

PARAMETER PI - 3.14159265 . TPI - 6.28318531 
PARAMETER XM - 2.5 , XNN - 3.0 
C 

C** Ci are solar activity variation variables 
C 

PARAMETER Cl - 383.0 . C2 - 3.32 . C3 - 1.80 
C 

C** Di are geomagnetic variation variables 
C 

PARAMETER DI - 28.0 , D2 - 0.03 . D3 - 1.0 , D4 - 100.0 . D5 - -0.08 
C 

C** Ei are semiannual variation variables 
C 

PARAMETER El - 2.41 . E2 - 0.349 , E3 - 0.206 , E4 - 6.2831853 
PARAMETER E5 - 3.9531708 . E6 - 12.5663706 . E7 - 4.3214352 
PARAMETER E8 - 0.1145 . E9 - 0.5 . E10 - 6.2831853 
PARAMETER Ell - 5.9742620 , El 2 - 2.16 

PARAMETER BETA - -0.6457718 , GAMMA - 0.7504916 . P - 0.1047198 
PARAMETER RE - 0.31 
C 

C** solar activity variation 
C 

TC - Cl + C2 * F10B + C3 • ( F10 - F10B ) 

C 

C** diurnal variation 
C 

ETA - 0.5 * ABS ( XLAT - SDA ) 

THETA - 0.5 * ABS ( XLAT + SDA ) 

TAU - SHA + BETA + P • SIN ( SHA + GAMMA ) 

IF ( TAU. GT. PI ) TAU - TAU - TPI 

IF ( TAU. LT.-PI ) TAU - TAU + TPI 

A1 - ( SIN ( THETA ) )**XM 

A2 - ( COS ( ETA ) )**XM 

A3 - ( COS ( TAU / 2. ) )**XNN 

B1 - 1.0 + RE • A1 

B2 - ( A2 — A1 ) / B1 

TV - B1 * ( 1 . + RE * B2 • A3 ) 

TL - TC * TV 

C 

C** geomagnetic variation 
C 

IF ( I1.EQ.1 ) THEN 

TG - DI • GI + D2 • EXP(GI) 

ELSE 

TG - D3 * GI + D4 * ( 1 - EXP ( D5 * GI ) ) 27 

END IF 
C 


o o o o o 


** semiannual variation 

G3 - 0.5 • ( 1.0 + SIN ( E10 • DY + Ell ) ) 

G3 - G3 ** El 2 

TAU1 - DY + E8 * ( G3 - E9 ) 

G1 - E2 + E3 • ( SIN ( E4 * TAU1 + E5 ) ) 

G2 - SIN ( E6 * TAU1+ E7 ) 

TS - El + F10B • G1 • G2 

** exospheric temperature 

TE - TL + TG + TS 

RETURN 

END 
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SUBROUTINE JAC ( Z , T . TZ . AN . A02 . AO , AA . AHE . AH , EM . 
DENS . DL ) 


••*•••**••**•******•*•••***•*• 

c** 

C** Subroutine ’JAC’ calculates the temperature TZ , the total density DENS 
C** and its logarithm DL, the mean molecular weight EM, the individual 
C** specie number densities for N, 02, 0, A, HE and H ( each preceded with 
C** an ’A’ ) at altitude Z given the exospheric temperature T. 

C** This subroutine uses the subroutine ’INTEGRATE’ and the function 
C** subprograms ’TEMP’ and *M0L_WT’. 

C** 

C** Rewritten by Mike Hickey, USRA 

C**************************************************************************** 

DIMENSION ALPHA(6) , El (6) . DI(6) , DIT(6) 

REAL* 4 M0L_WT 

PARAMETER AV - 6.02257E23 
PARAMETER QN - .78110 
PARAMETER 002 - .20955 
PARAMETER QA - .009343 
PARAMETER QHE - 1.289E-05 
PARAMETER RGAS - 8.31432 
PARAMETER PI - 3.14159265 
PARAMETER T0 - 183. 

GRAVITY ( ALTITUDE ) - 9.80665 / ( ( 1. + ALTITUDE / 6.356766E3 )**2 ) 


DATA ALPHA / 0.0 , 0.0 , 0.0 , 0.0 , -.380 . 0.0 / 

DATA El / 28.0134 . 31.9988 , 15.9994 , 39.948 , 4.0026 , 1.00797 / 


TX - 444.3807 + .02385 * T - 392.8292 * EXP ( -.0021357 * T ) 

A2 - 2. • (T-TX) / PI 

TX_T0 - TX - T0 

T1 - 1.9 • TX_T0 / 35. 

T3 - -1.7 • TX_T0 / ( 35.**3 ) 

T4 - -0.8 * TX_T0 / ( 35. * *4 ) 

TZ - TEMP ( Z , TX , T1 , T3 . T4 , A2 ) 


C** SECTION 1 
C** 

A - 90. 

D - AMIN1 ( Z . 105. ) 

FA - MOLJWT ( A ) * GRAVITY ( A ) / TEMP ( A , TX , T1 , T3 . T4 . A2 ) 
GD - GRAVITY ( D ) 

TD - TEMP ( D , TX , T1 , T3 , T4 , A2 ) 

EM - MOL_WT ( D ) 

FD - EM * GD / TD 

C Integrate gM/T from 90 to minimum of Z or 105 km :— 

CALL INTEGRATE ( A . D , FA . FD , R . TX , T1 . T3 . T4 . A2 ) 

C The number 2.1926E-8 » density x temperature/mean molecular weight at 90 km. 

DENS - 2. 1926E-8 • EM • EXP( -R / RGAS ) / TD 

FACTOR - AV • DENS 
PAR - FACTOR / EM 
FACTOR - FACTOR / 28.96 


C For altitudes below and at 105 km calculate the individual specie number 
C densities from the mean molecular weight and total density. 


IF ( Z. LE. 105 ) THEN 
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DL - A LOG 10 ( DENS ) 

AN - ALOG10 ( QN • FACTOR ) 

AA - ALOG10 ( QA • FACTOR ) 

AHE - ALOG10 ( QHE • FACTOR ) 

AO - ALOG10 ( 2. • PAR * ( 1 .-EM / 28.96 ) ) 

A02 - ALOG10 ( PAR • ( EM • ( 1.4002 ) / 28.96-1. ) ) 

AH - 0. 

C 

C*» Return to cal lino program 
C 

RETURN 


END IF 


C** SECTION 2 : This section is only performed for altitudes above 105 km 

C** 

C Note that having reached this section means that D in section 1 is 105 km. 

C 

C Calculate individual specie number densities from the total density and mean 
C molecular weight at 105 km altitude. 


DI(1) - QN • FACTOR 

DI(2) - PAR • (EM • (1.4002) / 28.96-1.) 
DI(3) - 2. • PAR • (1.- EM / 28.96) 

DI(4) - QA • FACTOR 
DI(5) - QHE • FACTOR 


FA1 - GD / TD 

FD1 - GRAVITY ( Z ) / TZ 

C Integrate g/T from 105 km to Z km :- 

CALL INTEGRATE ( D , Z , FA1 , FD1 , R . TX . T1 . T3 . T4 . A2 ) 


DO 41 1-1.5 

DIT(I) - DI (I) • ( TD / TZ ) **(1 .4ALPHA(I)) • EXP( -EI(I) • R / RGAS) 
IF ( DIT(I) . LE. 0. ) DIT(I) - 1.E-6 
41 CONTINUE 


C** This section calculates atomic hydrogen densities above 500 km altitude. 
C** Below this altitude , H densities are set to 10**— 6. 

C** SECTION 3 
C** 

IF ( Z .GT. 500. ) THEN 

A1 - 500. 

S - TEMP ( A1 , TX , T1 , T3 , T4 , A2 ) 

DI (6) - 10.** ( 73.13 - 39.4 • ALOG10 (S) + 5.5 • ALOG10(S) *ALOG10(S)) 

FaI - GRAVITY ( A1 ) / S 

CALL INTEGRATE ( A1 . Z . FAI , FD1 . R , TX . T1 , T3 . T4 . A2 ) 

DIT(6) - DI(6) • (S/TZ) • EXP ( -El (6) • R / RGAS ) 

ELSE 
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DIT (6) - 1.0 



END IF 


C For altitudes greater than 105 km , calculate total density and mean 
C molecular weight from individual specie number densities. 

DENS-0 

DO 42 1-1,6 

DENS - DENS + El (I) • DIT(I) / AV 

42 CONTINUE 

EM - DENS * AV / ( DIT(1 )+DIT(2)+DIT(3)+DIT(4)+DIT(5)+DIT(6) ) 
DL - ALOG10 (DENS) 

AN - ALOG10(DIT(1)) 

A02 - ALOG10(DIT(2)) 

AO - ALOG10(DIT(3)) 

AA - ALOG10(DIT(4)) 

AHE - ALOG10(DIT(5)) 

AH - ALOG1 0(DIT(6) ) 


RETURN 

END 
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FUNCTION TEMP ( ALT . TX , T1 , T3 . T4 , A2 ) 




C*« • * 
C** Function subprogram ’TEMP’ calculates the temperature at altitude ALT ** 
C** using equation (10) for altitudes between 90 and 125 km and equation ** 
C** (13) for altitudes greater than 125 km , from SAO Report 313. »* 
C*» ** 
C** Written by Mike Hickey. USRA ** 


PARAMETER BB - 4.5E-6 

U - ALT - 125. 

IF ( U .GT . 0. ) THEN 

TEMP - TX + A2 • ATAN ( T1 • U * ( 1 . + BB e (U**2.5)) / A2 ) 

ELSE 

TEMP - TX + T1 • U + T3 • (U**3) + T4 • (U**4) 

END IF 
END 


REAL FUNCTION M0L_WT*4 ( A ) 


***** 

c* • •* 

C** Subroutine ’MOU_WT* calculates the molecular weight for altitudes •• 

C*« between 90 and 105 km according to equation (1) of SAO report 313. ** 

C** Otherwise, MOL_WT is set to unity. ** 

C** *• 

C** Written by Mike Hickey. USRA ** 




DIMENSION B (7) 

DATA B / 28.15204 , -0.085586, 1.284E-4, -1.0056E-5. -1.021E-5, 
1.5044E-6. 9.9826E— 8 / 

IF ( A. GT. 105. ) THEN 

MO I WT - 1 . 

ELSE 

U - A - 100. 

MOL_WT - B (1) 

DO 1 1-2.7 

MO I WT - MOL_WT + B (I) • U ** ( 1-1 ) 

1 CONTINUE 

END IF 
END 
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SUBROUTINE INTEGRATE ( A . 0 . FA . FD , RR . TX . T1 , T3 , T4 , A2 ) 


C* 

C* 

C* 

C* 

C* 

C* 

C* 

C* 

C* 

C* 

c* 

c* 

c* 


Subroutine ’INTEGRATE’ performs the Simpson’s Rule quadrature ( SRQ4 ) 
implemented by G. F. Kuncir. 

A ■ lower limit of integration 

D - upper limit of integration 

EPS ■ relative error convergence criterion 

M - maximum number of integrations 

RR - result of integration 

N - number of integrations required to find R 
Written by Mike Hickey, USRA 


DIMENSION B (7) 

REAL*4 MOL_WT, R ( 0:10 ) 

GRAVITY ( ALTITUDE ) - 9.80665 / ( ( 1. + ALTITUDE / 6.356766E3 )**2 ) 
M - 10 


EPS - 

.0001 




SONE - 

( D - A 

) 

• 

( FA + FD ) * 0.5 

R (0) - 0.0 



DO 2 

N 

- 

1 , M 

NINT - 

2 •* N 




STOW - 

0. 




DEL - 

( D - A 

) 

/ 

FLOAT ( NINT ) 


DO 1 I - 1 , NINT . 2 
X - A + DEL • FLOAT ( I ) 

FX - MOL_WT (X) • GRAVITY ( X ) / TEMP ( X , TX . T1 , T3 . T4 , A2 ) 
STOW - STOW + FX 
1 CONTINUE 


R (N) - ( SONE + 4. * DEL • STOW ) / 3. 

C For altitudes greater than 180 km N usually reaches 5 or 6 before 
C convergence is achieved. Quite occasionally, R can appear to reach 
C convergence for N»3 at altitudes above 200 km. This anomalous and erroneous 
C result can be avoided by forcing N to go to values greater than 3 before 
C checking for convergence at altitudes greater than 180 km 

IF ( ( D .GE. 180. .AND. N .GE. 4 ).OR. D .LT. 180. ) THEN 

IF ( EPS • ABS ( R(N) ). GE. ABS ( R(N) - R(N-1) ) ) THEN 
RR - R (N) 

RETURN 

END IF 

END IF 

SONE - 0.25 • ( SONE + 3. • R (N) ) 

D IF ( N .EQ. M ) TYPE *. ’WARNING : Convergence not achieved’ 

2 CONTINUE 

RETURN 

END 
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SUBROUTINE SLV ( DEN . ALT , XLAT . DAY ) 

C******************«********************************************************* 

C»* Subroutine ’SLV’ computes the seasonal-latitudinal variation of density 
C** in the lower thermosphere in accordance with L. Jacchia, SAO 332, 1971. 
C** This affects the densities between 90 and 170 km. This subroutine need 
C** not be called for densities above 170 km, because no effect is observed. 
C** 

C** The variation should be computed after the calculation of density due to 
C** temperature variations and the density ( DEN ) must be in the form of a 
C** base 10 log. No adjustments are made to the temperature or constituent 
C** number densities in the region affected by this variation. 


C*e 



C** 

DEN 

■ densi ty ( 1 og10) 

C** 

ALT 

-alt! tude (km) 

C** 

XLAT 

• latitude (rad) 

C** 

DAY 

- day number 

c** 




C 


C** initialize density (DEN) ■ 0.0 
C 

DEN - 0.0 


C 

C** check if altitude exceeds 170 km 
C 

IF ( ALT. GT. 170. ) RETURN 
C 

C** compute density change in lower thermosphere 
C 


Z - ALT - 90. 

X - -0.0013 e Z • Z 
Y - 0.0172 • DAY + 1 .72 
P - SIN (Y) 

SP - ( SIN (XLAT) ) **2 
S - 0.014 • Z e EXP (X) 

D - S • P » SP 
C 

C** check to compute absolute value of 'XLAT' 
C 

IF ( XLAT. LT. 0. ) D - -0 
DEN - D 


RETURN 

END 
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SUBROUTINE SLVH ( DEN . DENHE . XLAT , SDA ) 




C** Subroutine ’SLVH' computes the seosono l-l ot i tud i nol variation of the 
C** helium number density according to L. Jacchia, SAO 332, 1971. This 
C** correction is not important below about 500 km. 


c»* 



c** 

DEN - 

density (Iog10) 

c** 

DENHE - 

helium number density (Iog10) 

c** 

XLAT - 

latitude (rad) 

c*» 

SDA - 

solar declination angle (rad) 


D0 - 10. ** DENHE 



A - ABS ( 0.65 * ( SDA / 0.40909079 ) ) 

B - 0.5 * XLAT 
C 

C** Check to compute absolute value of ’B’ 

C 

IF ( SDA. LT. 0. ) B - -B 
C 

C** compute X, Y, DHE and DENHE 
C 

X - 0.7854 - B 
Y - ( SIN (X) ) ** 3 
DHE- A • ( Y - 0.35356 ) 

DENHE - DENHE + DHE 
C 

C*« compute helium number density change 
C 

D1 - 10. ** DENHE 

DEL- D1 — D0 

RHO- 10. ** DEN 

DRHO - ( 6 . 646E-24 ) • DEL 

RHO - RHO + DRHO 

DEN - ALOG10 (RHO) 

RETURN 

END 
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SUBROUTINE FAIR5 ( DHEL1 ,DHEL2 ,DLG1 ,DLG2 ,IH , FDHEL , FDLG ) 


C»« This subroutine fairs between the region above 500 km, which invokes the 
C** seasono I — I at i tud i na I variation of the helium number density ( subroutine 
C** SLVH ), and the region below, which does not invoke any seasonal- 
C** latitudinal variation at all. 

C** 

C** INPUTS: DHEL1 

C** DHEL2 

C*» DLG1 

C** DLG2 

C*» IH 

C** IBFH 

C** OUTPUTS: FDHEL 
C** FDLG 

C** 

C** Written by Bill Jeffries, CSC, Huntsville, AL. 

C** ph. (205) 830-1000, x31 1 

C* ****************** ******************************* 


- helium number density before invoking SLVH 

- helium number density after invoking SLVH 

- total density before invoking SLVH 

- total density after invoking SLVH 
• height ( km ) — INTEGER 

- base fairing height ( km ) — INTEGER 

- faired helium number density 

- faired total density 


DIMENSION CZ ( 6 ) 

DATA CZ / 1.0, 0.9045085, 0.6545085, 0.3454915, 0.0954915, 0.0 / 
PARAMETER IBFH - 440 

C Height index 

I - ( IH - IBFH ) /10 + 1 
C Non-SLVH fairing coefficient 
CZI - CZ ( I ) 

C SLVH fairing coefficient 
SZI - 1.0 - CZI 
C Fai red densi ty 

FDLG - ( DLG1 • CZI ) + ( DLG2 • SZI ) 

C Faired helium number density 

FDHEL - ( DHEL1 • CZI ) + ( DHEL2 • SZI ) 

RETURN 

END 
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